%per bus kijken wat de spanning doet bij een uitval van een lijn.
%nu nog zonder vermogenvariatie
tic;
clear all;
IEEE_14_bus_original;
write_14bus;

aantal_lijnen = size(Line.con,1);
aantal_bussen = size(Bus.con,1);
% voltages_alle_bussen = zeros(aantal_lijnen+1,30);
% voltage_k_percent = zeros((aantal_lijnen+1)*9,30);
voltages_alle_bussen = zeros(aantal_lijnen,aantal_bussen);
%voltage_k_percent = zeros((aantal_lijnen)*7,30);
teller = 0;
for k = 1:1%0.7:0.05:1
    IEEE_14_bus_original;
%     PQ.con(:,4) = PQ.con(:,4)*k;
%     PQ.con(:,5) = PQ.con(:,5)*k;
    write_14bus;
    
    for jj=1:aantal_lijnen
        
        %Uitzetten netelement
        if jj~=1
            %terug inschakelen vorige netelement
            Line.con(jj-1,16)=1;
        end
        Line.con(jj,16)=0;
        write_14bus;
        %Berekenen Power Flow
            initpsat;
            Settings.freq = 60;
            Settings.absvalues= 'on';
            runpsat('IEEE_14_bus_modified','pf');
            voltages_alle_bussen(jj,:) = DAE.y(1+Bus.n:2*Bus.n)';
    end
%     Line.con(:,16)=1;
%         write_psatfile;
%         %Berekenen Power Flow
%             initpsat;
%             Settings.freq = 60;
%             Settings.absvalues= 'on';
%             runpsat('Nmin1','pf');
%             voltages_alle_bussen(aantal_lijnen+1,:) = DAE.y(1+Bus.n:2*Bus.n)';
            
    voltage_k_percent(aantal_lijnen*teller+1:aantal_lijnen*(teller+1),:) = voltages_alle_bussen; 
%     voltage_k_percent((aantal_lijnen+1)*teller+1:(aantal_lijnen+1)*(teller+1),:) = voltages_alle_bussen; 
    teller = teller+1;
    
end

%waar bevinden zich de dips?
dips_bussen_belasting = zeros(100,6);
teller = 0;

for l = 1:aantal_lijnen%((aantal_lijnen)*7)
    for ll = 1:aantal_bussen
        if voltage_k_percent(l,ll)<=0.9 && voltage_k_percent(l,ll)>=0.01
            teller = teller+1;
            %1e index van dips geeft geval weer, 2e index de bus
            dips_bussen_belasting(teller,1) = l;
            dips_bussen_belasting(teller,2) = ll;
            dips_bussen_belasting(teller,4) = voltage_k_percent(l,ll);
            
        end
    end
end
aantal_dips = teller;
%nullen verwijderen
for l=1:size(dips,1)
    if dips(l,1) == 0
        
        dips(l:size(dips,1), :) = [];
        break;
    end
end

%er zijn grenzen opgegeven van min en max voltage per belasting (0.9-1.1)
%checken of er een bus te laag of te hoog gaat:
lijst_getroffen_bussen = zeros(100,2);
teller = 1;
for l=1:size(dips)
    for k = 1:size(PQ.con)
    if dips(l,2)==PQ.con(k,1)
        lijst_getroffen_bussen(teller,2) = dips(l,2);
        lijst_getroffen_bussen(teller,1) = dips(l,1);
        teller = teller + 1;
    end
    end
end

%overzicht van meest gevoelige bussen
gevoeligheid_bussen = zeros(30,1);
for l = 1:size(dips,1)
    bus = dips(l,2);
    gevoeligheid_bussen(bus,1) = gevoeligheid_bussen(bus,1) + 1;
    
end

toc;


save lijn_uit dips_bussen_belasting aantal_dips 



